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Abstract. In this paper I review the theory and numerical simulations of non- 
linear dynamics of preheating, a stage of dynamical instability at the end of 
inflation during which the homogeneous infiaton explosively decays and deposits 
its energy into excitation of other matter fields. I focus on preheating in chaotic 
inflation models, which proceeds via broad parametric resonance. I describe 
a simple method to evaluate Floquet exponents, calculating stability diagrams 
of Mathieu and Lame equations describing development of instability in m^<j>'^ 
and preheating models. I discuss basic numerical methods and issues, and 
present simulation results highlighting non-equilibrium transitions, topological 
defect formation, late-time universality, turbulent scaling and approach to 
thermalization. I explain how preheating can generate large-scale primordial 
(non-Gaussian) curvature fluctuations manifest in cosmic microwave background 
anisotropy and large scale structure, and discuss potentially observable signatures 
of preheating. 
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1. Introduction 

The idea of inflation (a period of rapid quasi-exponential expansion of the Universe) 
neatly solves several long-standing issues in cosmology [TJ [2], and has been 
spectacularly confirmed by observations of the Cosmic Microwave Background (CMB) 
anisotropies [3j |4] . While the Universe is inflating, its contents is cold. But eventually, 
inflation has to end and the field driving the inflation must decay, depositing energy 
into high-energy particles. This process, known as reheating, "boils" the vacuum 
and starts the thermal history of the universe with the hot big bang. As universe 
continues to cool down, it could undergo more phase transitions, which would happen 
at symmetry breaking points of the theory. Very little is known about the fundamental 
physics at these energy scales, and cosmological observations could be our only source 
of information for the foreseeable future. No photons reach us directly from this 
epoch as the universe is filled with hot plasma and is opaque until recombination. 
Nevertheless, the expansion history of the early universe is imprinted on the sky in the 
form of primordial curvature fluctuations. With success of WMAP and the data from 
Planck soon to come, CMB observations are reaching precision required to disentangle 
other subdominant effects from Gaussian fluctuations due to simple inflation [5]. 

The most basic models of reheating involve inflaton decaying into one or more 
other scalar fields. Among the most interesting are the ones where decay is 
non-perturbative, for example proceeding through parametric resonance naturally 
happening in chaotic inflation models [51 [71 HI HI UHl [HI [H] , or tachyonic instability 
in hybrid inflation models [131 [Ml [TS] . For all their simplicity, these models 
have surprisingly rich physics involving non-equilibrium phase transitions. While 
initial stages of preheating are linear and instability development can be understood 
analytically [TTl [12] , dynamics could be chaotic [T6] , and the held evolution quickly 
becomes inhomogeneous and non-linear, so non-perturbative decay of the inflaton has 
to be studied numerically [ni[ll[li[20l[2Tl[221[2l[Ml[25l[2i[27]. In this paper I 
briefly review the theory of parametric resonance, go over general results of numerical 
simulations of non-linear field evolution during preheating, and discuss signatures of 
preheating that could potentially be observed. 

2. Analytical Theory of Preheating 

The basic model of reheating involves inflaton </) decaying into another scalar field x- 
The action describing two interacting scalar fields minimally coupled to gravity is 

s = y'|Y^-i5^''(</',M0,. + x.pX..)-^('/',x)}\/^d'x, (1) 

with potential 1^(0, x) containing the terms responsible for field masses and self- 
couplings, as well as their interaction. Polynomial field operators up to fourth order 
are renormalizable, so one usually takes ^ ni?(j3^ or i At/)'* infiaton potential for chaotic 
infiation, and \g'^4''^x^ coupling term [TT1[T2]. For simplicity, one can keep the decay 
field X massless. Couplings like \a4>x^ are also allowed and could be present [25] . 
although one would need a x^ self-interaction to keep the potential bounded from 
below. Models with various combinations of these potential bits have been studied in 
the literature; in this review, I will mainly focus on the one with quartic potential [12] 

v{c^,x) = \M' + \g^<l^\^- (2) 
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Figure 1. Dynamics of a single field chaotic inflation model with Xcj)"^ /i potential. 



In a flat homogeneous isotropic universe with Friedniann-Robertson- Walker metric 



ds^ = -dt^ +a{tfdy?, 
the field equations of motion are readily obtained from the action ([T|); they are 
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Hubble parameter H = a/a plays the role of friction term in field dynamics. Its value 
is determined by the (average) total energy density according to Friedmann equation 
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where the combined energy density of the two fields is 
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During chaotic inflation, potential energy of the inflaton causes Hubble friction to 
be large, and the motion of the fields is over-damped. The inflaton (j) slowly rolls 
down the potential until the damping becomes sub-critical, at which point it starts 
oscillating near the minimum of the potential with decreasing amplitude, as illustrated 
in Figure [T^. Evolution of the Hubble horizon size L = 1/H and the physical 
wavelength of comoving modes A = 27r a/fc is shown on Figure [TJ). During inflation L 
changes slowly, so that slow roll parameter e = ^''^^ ^ 1. When field is oscillating, 
Hubble horizon size grows as L cx according to the average equation of state, which 
is 1/3 for A^'* oscillator. Comoving modes stop exiting and begin re-entering the 
horizon when e = 1; this moment can be taken as the end of inflation. 
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Figure 2. Background oscillations of the inflaton a4> (black) and the zero mode 
of decay field ax (red) for values of coupling slightly inside (a) and outside (b) of 
the first resonance band in Figure [Sji. Amplitude of ax was scaled up by 3 ■ 10*. 



Other fields coupled to inflaton feel its oscillations through modulation of 
parameters in their equations of motion, such as the effective mass term g^cj)^ in 
equation ([5|). Periodic modulation can lead to parametric resonance and exponential 
growth of inhomogeneous excitations in the fields coupled to infiaton. This is a fairly 
generic feature of chaotic infiation models; let's see how this happens in our model ([2]). 

It is very useful to scale variables so that the inflaton oscillations are periodic and 
of constant amplitude. This is particularly easy in model (H]), which is conformally 
invariant apart from its coupling to gravity. Switching to conformal time drj = dt/a 
and scaling the field values according to their conformal weight, equation of motion for 
homogeneous infiaton (|H) becomes simply [acj))" + A(a0)'^ — a" cj) = 0. If one neglects 
the last term (which is small as a ~ ry), oscillating infiaton solution is 

'/'(^) = ^/(t), t = A^ $0(77-7/0) (8) 

where $0 is the amplitude of infiaton oscillations, and /(t) is a unit amplitude solution 
of the canonical anharmonic oscillator equation d"^ f /dr^+p = 0, which can be written 
exactly in terms of Jacobi elliptical cosine function, or its harmonic expansion |29| 

1 , 1 47r ^ cos(n — \)^t , , 

/(r)^cn(r,2-^) = 2^— ^ \ . (9) 

T cosh(n - 2)7r 

Function fir) is periodic with period T = 7r~^r^(-i), and its harmonic expansion 
is exponentially converging, so only a few terms are needed to accurately represent 
its shape. Substituting background infiaton solution (|8]) into equation of motion ([5]) 
for the coupled field x? and rescaling variables the same way we did for inflaton, one 
obtains the evolution equation for the Fourier mode of decay fleld Xk with comoving 
wavenumber k. In terms of rescaled parameters k = fc/(A2$o) and q = it is 

<P(.axk) 



+ 



(aXk) = 0. (10) 
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Figure 3. Stability diagram of Mathieu (a) and Lame (b) equations, using 
the same parametrization. White regions correspond to stable solutions, shaded 
regions are unstable, with brighter color corresponding to larger values of critical 
exponent ii (isolevels are spaced every A/x = 0.01185). 



Exact solution of the oscillating inflaton acj) is shown in Figure [21 along with 
homogeneous solution of the field ax coupled to it, for two slightly different values of 
the coupling .g^/A. As you can see, depending on the value of the coupling, evolution 
of the field x can be either exponentially unstable ([2ji), or merely oscillatory ([2b). 

Differential equations with periodic coefficients are studied in Floquet theory, 
and are often encountered in other branches of physics as well (for example, Bloch 
waves in condensed matter). Equation (|10l) . in particular, is known as Lame equation, 
while its counterpart for harmonic inflaton oscillations in m^cf)^ potential is Mathieu 
equation [3D]. According to Floquet 's theorem, equation ([TUj) admits solution of the 
form e'^^P(T), where /i is a complex number, and function P(t) is periodic. Floquet 
exponent fi depends on the parameters and g, and can be calculated by explicitly 
constructing such periodic function from the principal fundamental matrix solution 



W(r) 



Xiir) 
x'Ar) 



X2(t) 



(11) 



made up from two independent solutions xi and X2 with initial conditions W(0) = I. 
Integrating principal fundamental solution over a single period (which can be done 
efficiently and precisely using numerical integration) , and fixing coefficients in a linear 
combination ciXi(t) + 02X2(1") = e^^P(T) to satisfy P{T) = P(0) and P'{T) = P'(0), 
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Figure 4. Evolution of distributions of inflaton (left) and decay product (riglit) in 
a full non-linear simulation with coupling at the maximal resonance g^/A = 1.875. 



one finds that the value of Floquet exponent fi is given by 
e^'^ ^ QiT) + ^Q^{T)-W{T), 



(12) 



where Q and W are two invariants of the matrix W under a similarity transformation 
1 



Q 



■tr 



W = detW= 1. 



(13) 



The second invariant (Wronskian W) is conserved, so expression (IT^ simplifies to 



cosher = Q{T). 



(14) 



Stability diagrams of Mathieu and Lame equations showing contour plots of Re/i 
calculated in this fashion are presented on Figure [3l Wide bands in parameter space 
are unstable, with values of Floquet exponent ^ reaching as high as 0.237 for Lame 
equation. If the value of coupling lands you into one of these bands, inflaton will decay 
very efficiently into Xk particles for a wide range of wavenumbers k. This regime is 
known as a broad parametric resonance. 

Although not usually emphasized, the instability band structure of Lame and 
Mathieu equations is quite similar, as elliptic cosine © differs from cost by 4.3% 
total harmonic distortion and 18% larger period. The real difference between m^0^ 
and Xff)^ models is how parameters scale with expansion of the universe. While their 
values stay constant for A^"* model (which is nearly conformally invariant), for m?(tP' 



model expansion drags the values of parameters toward the stable point k 



0, 



eventually shutting off the resonance. Thus for preheating to be efficient in rn^fj)^ 
model, one needs a much larger initial value of q. 



3. Tackling Non-linear Evolution: Numerical Methods and Issues 

Broad parametric resonance amplifies quantum fluctuations of the fields, creating real 
particles in a state far from thermal equilibrium. Instability is exponentially rapid, 
and develops within a few dozen of inflaton oscillation (as illustrated in Figure H]), 
which is very fast on cosmological time scales. Once the energy density of created 
particles becomes comparable to that of the homogeneous inflaton, one can no longer 
treat the evolution perturbatively, and non-linearity of the coupling and back reaction 
of the created particles on the inflaton evolution has to be taken into account. The 
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Figure 5. Three-dimensional spatial discretization stencil (left) and summary of 
coefficients (rigfit) for minimal (top) and three isotropic discretization schemes. 



most straightforward way to do it is to solve field evolution equations numerically 

[HI HI]. 

Several codes are available for this purpose, most notably LATTICEEASY by 
Gary Felder and Igor Tkachev [JTl [32] which is widely used and modified, my own 
DEFROST [33] which offers improved performance and visualization capabilities, and 
a new GPU-accelerated CUDAEASY by Jani Sainio [34]. A pseudo-spectral code 
PSpectRE has just been released as well [35]. Most of the implementations (with the 
exception of PSpectRE) opt for a finite difference method to solve non-linear partial 
differential equations (|4|5p . The fields are discretized on a cubic spatial grid of spacing 
dx, and the spatial differential operators are approximated by finite differences 
D[X] ,^^,2 G[X] 



AX 



{vxy 



(15) 



(dx)'^ ' (dx)^ 
Discretizations of Laplacian operator involving only 26 nearest neighbours of a point 

x+l y+1 z+l 

^W=EEE^d(a)^a (16) 
1 1 z—1 



is going to be second order accurate, but truncation error can be made isotropic to 
fourth order 36 by taking coefficients Cq as summarized in Figure [5] A critical issue 
for accuracy of long-term cosmological simulations is discretization of gradient terms 
in energy density ([T]). Discretized energy is not necessarily conserved by discretized 
equations of motion, and gradient energy leaking off the grid affects equation of state 
and leads to large cumulative errors in expansion history of the universe [47l |49] . The 
best way to avoid this pitfall is to discretize Lagrangian ([T]) directly. One can show 
that the proper discretization of gradient square, variation of which leads to Laplacian 
discretization (|16p in equations of motion, is 

^ x+l y+1 z+l 

G[X]^-J2J2J2^^(-)(^-'^of- (17) 

x—1 y—1 z—1 



Once the spatial operators are discretized, time evolution problem becomes a system 
of coupled ordinary differential equations, which can be integrated using any of the 
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(a) t = 36.75 (b) t = 128.00 



Figure 6. Energy density p inside the simulation box soon after onset of 
instability (a) and during subsequent evolution (b) in preheating model ((2]|. 

usual methods. DEFROST uses leapfrog scheme, which is simple, fast, and second 
order accurate in time. Higher order scheme could be used if needed. Symplectic 
integrator developed in [JU] is capable of reaching machine precision levels, or hybrid 
integrator of Huang et. al. could be used if time operator splitting is not possible. 

4. Non-Linear Dynamics, Thermalization and Universality within Horizon 

Non-linear dynamics that soon takes over the evolution of scalar fields can be rather 
non-trivial. Preheating is essentially a non-equilibrium phase transition, and a lot 
of interesting things can happen. Over the years, detailed numerical studies have 
been carried out for many parametric resonance [171 Ull HOI HH [Ml HI] and tachyonic 
pn im [531 IM] preheating models. In this section, I highlight some features of non- 
linear dynamics in preheating on horizon scales. This is local physics as far as 
cosmology is concerned, as the Hubble horizon size at the end of chaotic inflation 
is tiny (roughly Im redshifted to the present day). 

Quantum fluctuations that fall into unstable bands of Figure [3] are amplified and 
become classical with large occupation numbers. Once linear instability develops, 
the field configuration becomes very inhomogeneous and particle description gets 
complicated by non-linear coupling. A useful tracer of field dynamics is the evolution of 
the total energy density ([T]), which is an adiabatic invariant for rapidly oscillating fields. 
Evolution of energy density p in preheating model ([2]) with g^/A = 1.875 is shown in 
Figure El For that value of the coupling, long wavelength fiuctuations grow fastest, 
leading to formation of large blobs as shown in Figure [6^. Once density contrast 
increases to about one, non-linear interaction kicks in and the density fragments to 
much smaller scales as shown in Figure [Sb- In particle description, this could be 
viewed as upscattering of modes to higher momenta by non-linear interaction term. 

Just as it happens in phase transitions, one can also produce topological defects 
during preheating, if the theory allows them. For example, in preheating model with 



Non-linear Dynamics and Primordial Curvature Perturbations from Preheating 9 




(a) t = 103.00 (b) t = 512.00 



Figure 7. Transient formation of topological defects (a) and late-time energy 
density configuration (b) in preheating model JTSj with global 0{2) symmetry. 



0(2)-syinmetric potential with a small field vacuum expectation value v 

V{ct>,x)^\K4>^+X^~v^f, (18) 

global cosmic strings can form [ini [20] • Initial instability in this model develops 
similarly to ([2]) with /\ = 2, but once the energy density dilutes enough to 
fill the ring at the bottom of potential, cosmic strings are produced. String cores 
(0^ + < w^/40) soon after formation are shown in Figure [7^. String loops within 
horizon are transient, and eventually will collapse and annihilate. After a while, energy 
density configuration once again reaches a highly fragmented state shown in Figure[7l3. 

An important open question is exactly how and when thermalization after 
preheating happens. Characteristic momentum of particles produced during linear 
stage of preheating is determined by the instability band structure. Typically, it is 
significantly less than that of thermal equilibrium value, so particles need to upscatter 
through non-linear interactions, and thermalization can be delayed by a long time 
[11] . This is supported by numerical simulations, which show slow scaling regime in 
evolution of field occupation numbers in model [37l |38] . 

While dynamics of the non-equilibrium phase transition can be rather varied, 
late stages of preheating appear to have a certain universality to them. After initial 
transient, the field evolution leads to a highly inhomogeneous state similar to that 
shown in FigureslHt) and[7}D which persists on long time scales, with slow fragmentation 
going on. A striking feature of this regime is that one-point probability distribution 
function of energy density contrast S = p/ p appears to be statistically stationary, and 
universal across a class of preheating models [33| . Figure [8] shows late-time energy 
density PDFs (with dilution due to expansion scaled out) for four different preheating 
models ending via broad parametric resonance. All of them fit lognormal distribution 

Pip)dp^'e.J-^^^^I^]'-^. (19) 




Figure 8. Universality of lognormal density distribution in various two-field 
preheating models with inflaton decaying via broad parametric resonance. 



It is very tempting to attribute scaling and universality of the late stages of preheating 
to scalar field turbulence [37l [38], especially since lognormal density distributions 
are known to arise in relativistic fluid turbulence |39) . but the subject needs further 
investigation. 



5. Large- Scale Primordial Fluctuations from Preheating 

As interesting as non-equilibrium field evolution could be, thermalization wipes out 
most of the details in the final state after the phase transition. However, dynamics of 
the transition could affect the expansion history of the universe, and leave an imprint in 
the observable large-scale curvature fluctuations produced during preheating [571 HH] ■ 
Inflation transforms sub-horizon quantum vacuum fluctuations in all the light 
fields into super-horizon classical fluctuations. These fluctuations are statistically 
homogeneous and isotropic Gaussian random fields, and are completely described 
by their spectra with amplitude P{k) ~ H'^/Att'^ evaluated at horizon crossing 
H = k/a. Causally disconnected patches on super- horizon scales evolve essentially 
independently, and large-scale curvature fluctuations $ are the difference SN in 
amount of expansion TV = In a different patches experience from the constant curvature 
hypersurface at the end of inflation to the constant density (and temperature) 
hypersurface once thermalization occurred [301 SIl HI] • Fluctuations of the inflaton 6(j) 
are usually the main source of metric curvature fluctuations with their amplitude 
enhanced by the slow- roll parameter P$(fc) — P^{k) / {2m^^t) . It is also entirely 
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Figure 9. Primordial fluctuation spectrum of field x produced by inflation |50| . 



possible to convert isocurvature modes from subdominant light fields into observable 
curvature perturbations, for example as it happens in curvaton-type scenarios |43[ 144] 
and modulated reheating |3B]. Resonant preheating dynamics can create and 
significantly amplify curvature fluctuations from isocurvature modes of light fields, as 
suggested by [171I1H] and calculated in [151 [50] , 

For simple preheating model ^ with small values of coupling g^/A, the second 
field X is light during inflation, and acquires fluctuation spectrum with power on 
super-horizon scales comparable to inflaton, as shown in Figure [5] Super-horizon 
fluctuations of x £^re converted to curvature fluctuations through preheating dynamics. 
The basic mechanism is that the flat x-direction of potential ([U is suddenly lifted due 
to expectation value of inhomogeneous terms like {54>^) when parametric resonance 
instability develops |49]. This will modulate equation of state based on the value of 
homogeneous mode in field x a-t a- time inhomogeneity develops, and create curvature 
fluctuations dependent on initial value of x on super-horizon scales. 

Calculating curvature fluctuations generated by preheating involves tracing 
minute differences in expansion history, from the end of inflation to thermalization, 
in non-linear three dimensional simulations of regions of the universe corresponding 
to different initial values of x on super-horizon scales. This is a very demanding 
numerical problem, not only in terms of computing power, but precision required. The 
first attempt encountered numerical difficulties [471 148) , and it required development 
of new numerical integration techniques to obtain the answer |49| . Skipping a lot of 
technical details which will be discussed elsewhere [SO] , the total curvature fluctuation 
$ produced in the inflation model ([2]) is 



where $g is the usual nearly Gaussian contribution from inflaton fluctuations (50, 
and the second uncorrelated term is generated by preheating from the super-horizon 




(20) 
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Figure 10. Non-linear transfer function -Fnl(x) connecting initial value of super- 
horizon mode Xini with curvature fluctuation &N it produces |49) . Thick red line 
shows result of averaging over substructure not resolved in CMB observations. 

mode of the field x- The exact distribution of the field x sampled on observable part 
of the sky depends on inflation history in extreme version of cosmic variance. The 
transfer function i^NL shown in Figure [10] is quite non-linear and could lead to non- 
Gaussian fluctuations of the form very different from the usual weak non-Gaussianity 
parametrization 

$(x) = $g(x) + /nl$gW. (21) 

The amplitude of curvature fluctuations produced by preheating in model ^ is 10~^, 
which is comparable to the curvature fluctuations from inflaton, so the two could 
potentially be disentangled by searching for non-Gaussian component in the observed 
CMB temperature anisotropy. 

6. Discussion: Going after Observable Signatures of Preheating 

Very little is known about how inflation actually ended, and what was the high 
energy physics like at those energy scales (or even what is the inflaton itself). 
Traces of reheating are hidden from us by opaque plasma in nearly thermal state, 
and are unobservable directly. One must seek signatures of preheating that survive 
thermalization and could be detected. These could include stable relics like topological 
defects [19, 20l[53j or primordial black holes [54j[55l|56], stochastic gravitational wave 
background produced by inhomogeneities during reheating [57l [58l [59l [60l EU [62] . 
or anomalies in the expansion history of the universe imprinted in the primordial 
curvature fluctuations [t5l[46ll47ll48ll49l[50l[5ll[52]. Of these, the last effect appears 
to be the most promising observationally, as stable relic formation is difficult without 
spoiling cosmology, and stochastic gravitational waves are very hard to detect. 

Curiously enough, the simple model of preheating ([2]) could generate non- 
Gaussian curvature fluctuations of observable amplitude which are intermittent, 
producing primordial "cold spots" \^ . Although arguable [1] , the CMB temperature 
map appears to have slight statistical anomalies in the form of a cold spot in 
the southern hemisphere |63j . and slight discrepancy in north-south temperature 
anisotropy spectra [64' . Exciting possibility is that primordial "cold spot" , whether 
from preheating or some other early universe source, could potentially explain both. 
A tell-tale signature of primordial non-Gaussian cold spot is the associated i?-mode 
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polarization pattern around it, which might be possible to test with Planck data j65) . 
Primordial potential "dips" would also manifest in formation of large scale structure. 
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